home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / nurb_polyg / NurbEval.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  6.3 KB  |  269 lines
  1. /*
  2.  * NurbEval.c - Code for evaluating NURB surfaces.
  3.  *
  4.  * John Peterson
  5.  */
  6.  
  7. #include <stdlib.h>
  8. #include <stdio.h>
  9. #include <math.h>
  10.  
  11. #include "nurbs.h"
  12. #include "drawing.h"
  13.  
  14.  
  15. /*
  16.  * Return the current knot the parameter u is less than or equal to.
  17.  * Find this "breakpoint" allows the evaluation routines to concentrate on
  18.  * only those control points actually effecting the curve around u.]
  19.  *
  20.  *    m   is the number of points on the curve (or surface direction)
  21.  *    k   is the order of the curve (or surface direction)
  22.  *    kv  is the knot vector ([0..m+k-1]) to find the break point in.
  23.  */
  24.  
  25. long
  26. FindBreakPoint( double u, double * kv, long m, long k )
  27. {
  28.     long i;
  29.  
  30.     if (u == kv[m+1])    /* Special case for closed interval */
  31.     return m;
  32.  
  33.     i = m + k;
  34.     while ((u < kv[i]) && (i > 0))
  35.     i--;
  36.     return( i );
  37. }
  38.  
  39. /*
  40.  * Compute Bi,k(u), for i = 0..k.
  41.  *  u        is the parameter of the spline to find the basis functions for
  42.  *  brkPoint    is the start of the knot interval ("segment")
  43.  *  kv        is the knot vector
  44.  *  k        is the order of the curve
  45.  *  bvals    is the array of returned basis values.
  46.  *
  47.  * (From Bartels, Beatty & Barsky, p.387)
  48.  */
  49.  
  50. static void
  51. BasisFunctions( double u, long brkPoint, double * kv, long k, double * bvals )
  52. {
  53.     register long r, s, i;
  54.     register double omega;
  55.  
  56.     bvals[0] = 1.0;
  57.     for (r = 2; r <= k; r++)
  58.     {
  59.     i = brkPoint - r + 1;
  60.     bvals[r - 1] = 0.0;
  61.     for (s = r-2; s >= 0; s--)
  62.     {
  63.         i++;
  64.         if (i < 0)
  65.         omega = 0;
  66.         else
  67.         omega = (u - kv[i]) / (kv[i + r - 1] - kv[i]);
  68.         bvals[s + 1] = bvals[s + 1] + (1 - omega) * bvals[s];
  69.         bvals[s] = omega * bvals[s];
  70.     }
  71.     }
  72. }
  73.  
  74. /*
  75.  * Compute derivatives of the basis functions Bi,k(u)'
  76.  */
  77. static void
  78. BasisDerivatives( double u, long brkPoint, double * kv, long k, double * dvals )
  79. {
  80.     register long s, i;
  81.     register double omega, knotScale;
  82.  
  83.     BasisFunctions( u, brkPoint, kv, k - 1, dvals );
  84.  
  85.     dvals[k-1] = 0.0;        /* BasisFunctions misses this */
  86.  
  87.     knotScale = kv[brkPoint + 1L] - kv[brkPoint];
  88.  
  89.     i = brkPoint - k + 1L;
  90.     for (s = k - 2L; s >= 0L; s--)
  91.     {
  92.     i++;
  93.     omega = knotScale * ((double)(k-1L)) / (kv[i+k-1L] - kv[i]);
  94.     dvals[s + 1L] += -omega * dvals[s];
  95.     dvals[s] *= omega;
  96.     }
  97. }
  98.  
  99. /*
  100.  * Calculate a point p on NurbSurface n at a specific u, v using the tensor product.
  101.  * If utan and vtan are not nil, compute the u and v tangents as well.
  102.  *
  103.  * Note the valid parameter range for u and v is
  104.  * (kvU[orderU] <= u < kvU[numU), (kvV[orderV] <= v < kvV[numV])
  105.  */
  106.  
  107. void
  108. CalcPoint(double u, double v, NurbSurface * n,
  109.        Point3 * p, Point3 * utan, Point3 * vtan)
  110. {
  111.     register long i, j, ri, rj;
  112.     register Point4 * cp;
  113.     register double tmp;
  114.     register double wsqrdiv;
  115.     long ubrkPoint, ufirst;
  116.     double bu[MAXORDER], buprime[MAXORDER];
  117.     long vbrkPoint, vfirst;
  118.     double bv[MAXORDER], bvprime[MAXORDER];
  119.     Point4 r, rutan, rvtan;
  120.  
  121.     r.x = 0.0;
  122.     r.y = 0.0;
  123.     r.z = 0.0;
  124.     r.w = 0.0;
  125.  
  126.     rutan = r;
  127.     rvtan = r;
  128.  
  129.     /* Evaluate non-uniform basis functions (and derivatives) */
  130.  
  131.     ubrkPoint = FindBreakPoint( u, n->kvU, n->numU-1, n->orderU );
  132.     ufirst = ubrkPoint - n->orderU + 1;
  133.     BasisFunctions( u, ubrkPoint, n->kvU, n->orderU, bu );
  134.     if (utan)
  135.     BasisDerivatives( u, ubrkPoint, n->kvU, n->orderU, buprime );
  136.  
  137.     vbrkPoint = FindBreakPoint( v, n->kvV, n->numV-1, n->orderV );
  138.     vfirst = vbrkPoint - n->orderV + 1;
  139.     BasisFunctions( v, vbrkPoint, n->kvV, n->orderV, bv );
  140.     if (vtan)
  141.     BasisDerivatives( v, vbrkPoint, n->kvV, n->orderV, bvprime );
  142.  
  143.     /* Weight control points against the basis functions */
  144.  
  145.     for (i = 0; i < n->orderV; i++)
  146.     for (j = 0; j < n->orderU; j++)
  147.     {
  148.         ri = n->orderV - 1L - i;
  149.         rj = n->orderU - 1L - j;
  150.  
  151.         tmp = bu[rj] * bv[ri];
  152.         cp = &( n->points[i+vfirst][j+ufirst] );
  153.         r.x += cp->x * tmp;
  154.         r.y += cp->y * tmp;
  155.         r.z += cp->z * tmp;
  156.         r.w += cp->w * tmp;
  157.  
  158.         if (utan)
  159.         {
  160.         tmp = buprime[rj] * bv[ri];
  161.         rutan.x += cp->x * tmp;
  162.         rutan.y += cp->y * tmp;
  163.         rutan.z += cp->z * tmp;
  164.         rutan.w += cp->w * tmp;
  165.         }
  166.         if (vtan)
  167.         {
  168.         tmp = bu[rj] * bvprime[ri];
  169.         rvtan.x += cp->x * tmp;
  170.         rvtan.y += cp->y * tmp;
  171.         rvtan.z += cp->z * tmp;
  172.         rvtan.w += cp->w * tmp;
  173.         }
  174.     }
  175.  
  176.     /* Project tangents, using the quotient rule for differentiation */
  177.  
  178.     wsqrdiv = 1.0 / (r.w * r.w);
  179.     if (utan)
  180.     {
  181.     utan->x = (r.w * rutan.x - rutan.w * r.x) * wsqrdiv;
  182.     utan->y = (r.w * rutan.y - rutan.w * r.y) * wsqrdiv;
  183.     utan->z = (r.w * rutan.z - rutan.w * r.z) * wsqrdiv;
  184.     }
  185.     if (vtan)
  186.     {
  187.     vtan->x = (r.w * rvtan.x - rvtan.w * r.x) * wsqrdiv;
  188.     vtan->y = (r.w * rvtan.y - rvtan.w * r.y) * wsqrdiv;
  189.     vtan->z = (r.w * rvtan.z - rvtan.w * r.z) * wsqrdiv;
  190.     }
  191.  
  192.     p->x = r.x / r.w;
  193.     p->y = r.y / r.w;
  194.     p->z = r.z / r.w;
  195. }
  196.  
  197. /*
  198.  * Draw a mesh of points by evaluating the surface at evenly spaced
  199.  * points.
  200.  */
  201. void
  202. DrawEvaluation( NurbSurface * n )
  203. {
  204.     Point3 p, utan, vtan;
  205.     register long i, j;
  206.     register double u, v, d;
  207.     SurfSample ** pts, *sptr;
  208.  
  209.     long Granularity = 10;  /* Controls the number of steps in u and v */
  210.  
  211.     /* Allocate storage for the grid of points generated */
  212.  
  213.     CHECK( pts = (SurfSample**) malloc( (Granularity+1L) * sizeof( SurfSample* )));
  214.     CHECK( pts[0] = (SurfSample*) malloc( (Granularity+1L)*(Granularity+1L)
  215.             * sizeof( SurfSample )));
  216.     for (i = 1; i <= Granularity; i++)
  217.     pts[i] = &(pts[0][(Granularity+1L) * i]);
  218.     sptr = pts[0];
  219.  
  220.     /* Compute points on curve */
  221.  
  222.     for (i = 0; i <= Granularity; i++)
  223.     {
  224.     v = ((double) i / (double) Granularity)
  225.         * (n->kvV[n->numV] - n->kvV[n->orderV-1])
  226.         + n->kvV[n->orderV-1];
  227.  
  228.     for (j = 0; j <= Granularity; j++)
  229.     {
  230.         u = ((double) j / (double) Granularity)
  231.         * (n->kvU[n->numU] - n->kvU[n->orderU-1])
  232.         + n->kvU[n->orderU-1];
  233.  
  234.         CalcPoint( u, v, n, &(pts[i][j].point), &utan, &vtan );
  235.         (void) V3Cross( &utan, &vtan, &p );
  236.         d = V3Length( &p );
  237.         if (d != 0.0)
  238.         {
  239.         p.x /= d;
  240.         p.y /= d;
  241.         p.z /= d;
  242.         }
  243.         else
  244.         {
  245.         p.x = 0;
  246.         p.y = 0;
  247.         p.z = 0;
  248.         }
  249.         pts[i][j].normLen = d;
  250.         pts[i][j].normal = p;
  251.         pts[i][j].u = u;
  252.         pts[i][j].v = v;
  253.     }
  254.     }
  255.  
  256.     /* Draw the grid */
  257.  
  258.     for (i = 0; i < Granularity; i++)
  259.     for (j = 0; j < Granularity; j++)
  260.     {
  261.         (*DrawTriangle)( &pts[i][j], &pts[i+1][j+1], &pts[i+1][j] );
  262.         (*DrawTriangle)( &pts[i][j], &pts[i][j+1],     &pts[i+1][j+1] );
  263.     }
  264.  
  265.     free( pts[0] );
  266.     free( pts );
  267. }
  268.  
  269.